function [ys,params,check] = private_1s2c_perfect2_varyTaux_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% computes the steady state for the NK_baseline.mod and uses a numerical

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;

%%%%%%%%%% Don't change the above code %%%%%%%%%%%%
global theta dbar  L gL alpha rbest rho sigmaH sigmaF ns n betas d alphas fspace_smolyak ilinear ilog_rL delta beta


% private eqm
r_private = rbest;
n  = 2;
ns = 1;
ilog_rL = 1; % use exp(rL)/(1+exp(rL))

dij = dbar;
d=dij*ones(n); %dij is (1+iceberg) cost producing in j selling to i
d=d-(dij-1)*eye(n);


for i=1:n
    eval(['L(',num2str(i),',1)=L',num2str(i),';'])
    eval(['alpha(',num2str(i),')=alpha',num2str(i),';'])
end


%% solve for the steady state of gov assuming dG2/dT=0!!!
% This is not the true steady state, it's only for us to get initial guess
% --- steady state ---

options = optimoptions('fsolve', 'TolFun', 1e-12);
options.Display = 'off';

xguess(1:2) = 1.1;
xguess(3) = 0.2;

[xval,fval_sol,eflag_ss] = fsolve(@(xval) solve_1s_2c_gov_ss_nogL(xval),xguess,options); %iter


[Tss,xss,Pss,pi_ss, fval_ss] = eval_1s_2c_gov_ss_nogL(xval);

xval_ss= xval;
wss = xval(1:n);
extax_ss = xval(3);
taf_ss = 0;
rss = rbest;

taux = extax_ss;


uc = xss(1)^(-sigmaH);

               
gammaT_ss = wss(2)/alpha(2)*(1+theta)/theta*pi_ss(2,1)*uc*extax_ss/(1+extax_ss);
mu = gammaT_ss*alpha(2);

%% write to our variable names in dynare
for j=1:ns
    eval(['wedge',num2str(j),'=0;'])
end

rL = rbest*L;
Lp_ss=(1-rbest)*L;


eval(['bterm= 0;'])
for i=1:n
    eval(['x',num2str(i),'=log(xss(',num2str(i),'));'])
    eval(['w',num2str(i),'=log(wss(',num2str(i),'));'])
    eval(['P',num2str(i),'=log(Pss(',num2str(i),'));'])
    eval(['riskfree',num2str(i),'=log(1/beta);']);
    if (i>1)
        eval(['mu',num2str(i),'=mu;'])
    end


    eval(['T',num2str(i),'=log(Tss(',num2str(i),'));'])
    eval(['r',num2str(i),'=log(rbest);'])
    eval(['Lr',num2str(i),'=log(rL(',num2str(i),'));'])
    eval(['Lp',num2str(i),'=log(Lp_ss(',num2str(i),'));'])

    if (i>1)

        eval(['taux',num2str(i),'= extax_ss;'])
        eval(['taf',num2str(i),'= taf_ss;'])


    end


    for ii=1:n
        eval(['pi',num2str(i),num2str(ii),'= pi_ss(',num2str(i),',',num2str(ii),');'])
    end



end

%% end own model equations, below from dynare

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
